function eq = solveEquilibriumDebtOverhang(z0,lsap,gamma,b)
% Solve the equilibrium for the section with debt overhang and firm insolvency
% This is similar to the baseline case with the difference that
% Productivity is now endogenous and given by z0 \bar{S}
% where \bar{S} =gamma + (1-gamma) *F(z0p/(1-exp(-rho)))
% Here: gamma captures the efficiency of the bankruptcy system
% and F(.) captures the distribution of initial debt/cash across firms
% 
% We assume F(.) is uniformly distributed over [-b,b]
% Therefore, F(\bar{b}) = (\bar{b}+b)/2b

% Inputs:
% lsap = etag*(wg-1); % Size of the lsap program
% gamma = Efficiency of the bankruptcy system
% b = Highest outstanding debt level: Assume uniform distribution over [-b,b]

% Returns: An equilibrium structure with (among other things):
% The equilibrium price p and the interest rate rf (if multiple equilibria, return the best equilibrium)
% The required Sharpe ratio curve (for plotting purposes)
% The actual Sharpe ratio curve with rf=0 (for plotting purposes)

assert(lsap>=0);  % We are only interested in cases with positive lsap
assert(lsap<1); % The case above 1 doesn't make much sense
assert(z0<=1);  % We are only interested in cases with productivity shocks are below baseline
assert(gamma>=0 & gamma<=1); 
assert(b>=0);
eq.lsap = lsap;
eq.z0 = z0;
eq.gamma =gamma;
eq.b = b;

% Underlying parameters(loaded through globals)
global g rho sigma PStar taub tauh k l tau tauBarBench

ph = exp(rho + g - sigma^2*(1-lsap)/tauh);
tau = taub*k+tauh*(1-k);
% Store all parameters for later access
eq.params.g = g; eq.params.rho = rho; eq.params.sigma = sigma;
eq.params.tauh =tauh; eq.params.taub = taub;
eq.params.tau = tau;

eq.ph=ph;
% Baseline price without a demand recession
% Will adjust this for cases with demand recession
p = 1;

if ph >= 1
    % Simple case where households are sufficiently risk tolerant (with lsap)
    % No demand recession (regardless of z0 or bankruptcies)
    eq.unique=1;
    p=1;
    % No need to change the baseline levels
else
% ph < 1
% There can be a demand recession

% Focus on the case in which there is a unique interior equilibrium
% (Change the parameters appropriately to ensure this)

% Set the appropriate arrays (for plotting purposes)
pArr = linspace(ph,1,100); % Normalized price array
eq.pArr = pArr;

sharpeActualArr = (g + rho - log(pArr))/sigma;
eq.sharpeActualArr = sharpeActualArr;

% Calculate banks' implied balance sheets, and the implied risk tolerance, for each p
SArr = min(1,(z0*pArr/(1-exp(-rho))+b)/(2*b)); % Fraction of solvent firms 
% For the baseline, b=0, the right side gives Inf and this reduces to Sarr=1 (no insolvencies, which is correct)

SBarArr = gamma + (1-gamma)*SArr; % Effective productivity
alphaArr = max(0,k*(1-l./(SBarArr.*z0.*pArr)));
eq.alphaArr = alphaArr;
tauBarArr = alphaArr*taub+(1-alphaArr)*tauh;
eq.tauBarArr = tauBarArr;

sharpeRequiredArr = sigma*(1-lsap)./tauBarArr;
eq.sharpeRequiredArr = sharpeRequiredArr;

% Find the unique intersection (assuming parameters are such that the
% equilibrium is unique)
   [sol,numCross] = findCross(sharpeRequiredArr,sharpeActualArr,pArr);
   assert(numCross==1); % Intersection is supposed to be unique
   p = sol{1}.x;
   ind = sol{1}.xind;
   % eq.Sharpe = sol.y;
  
end

% Set the key equilibrium output
% The key thing that we should get right from the earlier analysis is p
% Given p (if correctly calculated), we can recover all else
S = min(1,(z0*p/(1+exp(-rho))+b)/(2*b)); % Fraction of solvent firms 
SBar = gamma + (1-gamma)*S; % Effective productivity

alpha = max(0,k*(1-l/(SBar*z0*p)));
tauBar = alpha*taub+(1-alpha)*tauh;
rf = max(0,rho + g - (1-lsap)*sigma^2/tauBar); % Accommodates corner solutions
eq.p = p;
eq.ind = ind;
eq.S = S;
eq.SArr = SArr;
eq.SBar = SBar;
eq.SBarArr = SBarArr;
eq.alpha= alpha;
eq.tauBar = tauBar;
eq.rf = rf;


function [sol,numCross] = findCross(y1Arr,y2Arr,xArr)
    % To accommodate corner solutions, choose y1Arr and y2Arr so that
    % if y1Arr<y2Arr throughout then xArr(end) is the solution (and vice versa)
    
    % Find all crossings of two relations
    % y1(x1) and y2(x2)
    % defined on the same array, xArr
    % Requirement: The relations are continuous so every crossing will be interpreted as an equilibrium
    % Requirement: y1Arr(xArr)
    % Output:
    % sol{i}.y and sol{i}.x have the solutions corresponding to crossing i
    % sol{i}.xind has the x-index corresponding to the solution with crossing i
    % numCross: Number of equilibria      
    %
    % Corner cases: If there is no crossing, assume a corner solution (depending on y1<y2 or y1>y2) 
    numCross =0; 
    ind1Arr = []; 
    lastDif = 0; % Keeps track of the last difference, y1-y2, at the last step
    % newDif will keep track of the difference at this iteration
for ind=1:length(xArr)
    x = xArr(ind);
    y1 = y1Arr(ind);
    y2 = interp1(xArr,y2Arr,x,'linear');
    newDif = y1-y2;
    if isnan(y2)
        % No intersection
        continue;
    elseif newDif==0
        display('Crossing found!');
        numCross = numCross+1;
        sol{numCross}.y = y1;
        sol{numCross}.x = x;
        sol{numCross}.xind = ind;
    elseif newDif<0 & lastDif>0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        % Find the exact solution
        numCross = numCross+1;
        % Do a linear interpolation based on the magnitude of the differences
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    elseif newDif>0 & lastDif<0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        numCross = numCross+1;
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    end
    lastDif =newDif;
end  

    % If at the end of the loop numEq is zero, make some adjustment for corner solution
    if numCross ==0
        if lastDif<0
            % This means that y1<y2 throughout
            % Set the highest level of xArr as a corner solution 
            numCross=1;
            sol{numCross}.y = y1Arr(end);
            sol{numCross}.x = xArr(end);
            sol{numCross}.xind = length(xArr);
        elseif lastDif >0
            % This means that y1>y2 throughout
            % Set the lowest level of xArr as a corner solution 
            numCross= 1;
            sol{numCross}.y = y1Arr(1);
            sol{numCross}.x = xArr(1);
            sol{numCross}.xind = 1;
        end
    end

end

end